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ABSTRACT 

We analyze 3D numerical simulations of driven incompressible magnetohydro- 
dynamic (MHD) turbulence in a periodic box threaded by a moderately strong 
external magnetic field. We sum over nonlinear interactions within Fourier wave- 
bands and find that the time scale for the energy cascade is consistent with the 
Goldreich-Sridhar model of strong MHD turbulence. Using higher order lon- 
gitudinal structure functions we show that the turbulent motions in the plane 
perpendicular to the local mean magnetic field are similar to ordinary hydrody- 
namic turbulence while motions parallel to the field are consistent with a scaling 
correction which arises from the eddy anisotropy. We present the structure tensor 
describing velocity statistics of Alfvenic and pseudo-Alfvenic turbulence. Finally, 
we confirm that an imbalance of energy moving up and down magnetic field lines 
leads to a slow decay of turbulent motions and speculate that this imbalance is 
common in the interstellar medium where injection of energy is intermittent both 
in time and space. 

Subject headings: ISM:general-MHD-turbulence 

1. Introduction 

The interstellar medium (ISM) is complicated and dynamic. The magnetic field and 
dynamic pressure (pv 2 /2) usually dominate the thermal pressure (nkT), dramatically influ- 
encing the star formation rate (see McKee 1999 for a review). There are cosmic rays that 
provide pressure and heating as well. 
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One approach to studying the ISM is to perform time-dependent numerical simulations 
to model the ISM, including as many of the interacting phenomena as practical. Of course, 
the physics included in such models must necessarily be highly simplified, and it is difficult 
to determine which features of the final model result from which physical assumptions (or 
initial conditions). Our approach is to use simplified numerical simulations to study the 
influences of various physical phenomena in isolation. We want to obtain a physical feeling 
for the general effects that each phenomenon has on the nature of the ISM. In this paper 
we will consider the influence of random forces per unit volume on MHD turbulence in 
an incompressible medium. It is obvious that the real ISM is compressible, but we want to 
separate the effects of magnetic turbulence from those involving compression. In later papers 
we will include compression for comparison with the present models, thereby isolating its 
importance directly. 

Historically hydrodynamic turbulence in an incompressible fluid was successfully de- 
scribed by the eddy cascade (Kolmogorov 1941), but MHD turbulence was first modeled 
by wave turbulence (Iroshnikov 1963, Kraichnan 1965; hereinafter IK). This theory assumes 
isotropy of the energy cascade in Fourier space, an assumption which has attracted severe 
criticism (Montgomery & Turner 1981; Shebalin et al. 1983; Montgomery & Matthaeus 1995; 
Sridhar & Goldreich 1994). Indeed, the magnetic field defines a local symmetry axis since it 
is easy to mix field lines in directions perpendicular to the local B and much more difficult 
to bend them. The idea of an anisotropic (perpendicular) cascade has been incorporated 
into the framework of the reduced MHD approximation (Strauss 1976; Rosenbluth 1976; 
Montgomery 1982; Zank & Matthaeus 1992; Bhattacharjee, Ng, & Spangler 1998). 

In a turbulent medium, the kinetic energy associated with large scale motions is greater 
than that of small scales. However, the strength of the local mean magnetic field is almost 
the same on all scales. Therefore, it becomes relatively difficult to bend magnetic field 
lines as we consider smaller scales, leading to more pronounced anisotropy. A self-consistent 
model of MHD turbulence which incorporates this concept of scale dependent anisotropy 
was introduced by Goldreich & Sridhar (1995) (henceforth GS95). 

Within the GS95 theory the energy cascade becomes anisotropic as a consequence of 
the resonant conditions for 3-wave interactions. A strict application of the resonant 3-wave 
interaction conditions gives an energy cascade which is purely in the direction perpendicular 
to the external field. However, it is intuitively clear that the increase in k± must at some 
point start affecting k\\. 

The cornerstone of the GS95 theory is the concept of a 'critically balanced' cascade, 
where h\\V~A ~ k±vi, where k± and k\\ are wave numbers perpendicular and parallel to the 
background field an v\ is the r.m.s. speed of turbulence at the scale I. In this model, the 
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Alfven rate (A^Va) is equal to the eddy turnover rate (k±vi). Using this concept, Goldreich 
& Sridhar showed that the energy cascade is not strictly perpendicular to the background 
field, but is relaxed so that k» oc k 2 / 3 . Their model predicts that the one-dimensional 
energy spectrum is Kolmogorov-type if expressed in terms of perpendicular wavenumbers, 
i.e. E(k±) oc k^ 3 . 

Numerical simulations by Cho & Vishniac (2000a, henceforth CV00) and Maron & 
Goldreich (2001, hereinafter MG01) have mostly supported the GS95 model and helped to 
extend it. Both analyses stressed the point that scale dependent anisotropy can be measured 
only in local coordinate frames which are aligned with the locally averaged magnetic field 
direction. Cho & Vishniac (2000a) calculated the structure functions of velocity and magnetic 
field in the local frames, and found that the contours of the structure functions do show 
scale dependent anisotropy, consistent with the predictions of the GS95 model. In their 
calculation, the strength of the uniform background magnetic field is roughly the same as the 
the r.m.s. velocity. MG01 tested the GS95 model for a much stronger uniform background 
field and also obtained results supporting the GS95 model, but they produced E(K) oc k~ 3 / 2 . 
They also calculated time scales of turbulence, interactions between pseudo- and shear- 
Alfvenic modes, growth of imbalance, and intermittency. Other related recent numerical 
simulations include Matthaeus et al. (1998), Miiller & Biskamp (2000), and Milano et al. 
(2001). 

These studies left a number of unresolved issues, including the exact scaling relations, 
the comparison of intermittency in MHD and in hydrodynamic turbulence, and the time 
scale of turbulence decay. Moreover, for many practical applications a more quantitative 
description of MHD turbulence statistics is necessary. These are vital for understanding var- 
ious astrophysical processes, including star formation (McKee 1999), cosmic ray propagation 
(Kota & Jokipii 2000), and magnetic reconnection (Lazarian & Vishniac 1999). 

In this paper, we further investigate implications of the GS95 model. In §2, we explain 
our numerical method. In §3, we further elucidate the scaling relation implied by the GS95 
model. In particular, we discuss the time scale, velocity scaling relations, and intermittency. 
In §4, we derive the correlation tensor and discuss some astrophysical applications. While 
the GS95 model predicts the MHD turbulence decays in just one eddy turnover time, in §5, 
we show that the decay time scale increases when the cascade is unbalanced and discuss 
some consequences of this fact. In §6, we briefly discuss the implications of this work. In §7 
we give a summary and our conclusions. As before, we consider the case where the uniform 
background magnetic field energy density is comparable to the turbulent energy density. 
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2. Method 

2.1. Numerical Method 

We have calculated the time evolution of incompressible magnetic turbulence subject to 
a random driving force per unit mass. We have adopted a pseudospectral code to solve the 
incompressible MHD equations in a periodic box of size 2ir: 

— = (V x v) x v - (V x B) x B + z/V 2 v + f + VP', (1) 

- = Vx(vxB) + J) V 2 B, (2) 

V • v = V • B = 0, (3) 

where f is a random driving force, P' = P/p + v • v/2, v is the velocity, and B is magnetic 
field divided by (Anp) 1 ^ 2 . In this representation, v can be viewed as the velocity measured 
in units of the r.m.s. velocity, v, of the system and B as the Alfven speed in the same units. 
The time t is in units of the large eddy turnover time (~ L/v) and the length in units of L, 
the inverse wavenumber of the fundamental box mode. In this system of units, the viscosity 
v and magnetic diffusivity rj are the inverse of the kinetic and magnetic Reynolds numbers 
respectively. The magnetic field consists of the uniform background field and a fluctuating 
field: B = B + b. We use 21 forcing components with 2 < k < \^12, where wavenumber k 
is in units of L -1 . Each forcing component has correlation time of one. The peak of energy 
injection occurs at k ps 2.5. The amplitudes of the forcing components are tuned to ensure 
v « 1 We use exactly the same forcing terms for all simulations. The Alfven velocity of the 
uniform background field, B Q , is set to 1. We consider only cases where viscosity is equal to 
magnetic diffusivity: 

v = i). (4) 

In pseudo spectral methods, the temporal evolution of equations (1) and (2) are followed in 
Fourier space. To obtain the Fourier components of nonlinear terms, we first calculate them 
in real space, and transform back into Fourier space. The average kinetic helicity in these 
simulations is not zero. However, previous tests have shown that our results are insensitive 
to the value of the kinetic helicity. In incompressible fluid, P' is not an independent variable. 
We use an appropriate projection operator to calculate VP' term in Fourier space and also 
to enforce divergence- free condition (V • v = V ■ B = 0). We use up to 256 3 collocation 
points. We use an integration factor technique for kinetic and magnetic dissipation terms 
and a leap-frog method for nonlinear terms. We eliminate the 2 At oscillation of the leap-frog 
method by using an appropriate average. At t — 0, the magnetic field has only its uniform 
component and the velocity field is restricted to the range 2 < k < 4 in wavevector space. 
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Hyperviscosity and hyperdiffusivity are used for the dissipation terms (see Table 1). 
The power of hyperviscosity is set to 8, so that the dissipation term in the above equation 
is replaced with 

-^ 8 (V 2 ) 8 v, (5) 

where v% is determined from the condition Uh(N/2) 2h At ~ 0.5 (see Borue and Orszag 1996) 
Here At is the time step and N is the number of grid points in each direction. The same 
expression is used for the magnetic dissipation term. We list parameters used for the simu- 
lations in Table 1. We use the notation XY-B Z, where X = 256, 144 refers to the number 
of grid points in each spatial direction; Y = H refers to hyperviscosity; Z=l refers to the 
strength of the external magnetic field. 

Diagnostics for our code can be found in Cho and Vishniac (2000b). For example, our 
code conserves total energy very well in simulations with v = i] = and the average energy 
input (=f • v) is almost exactly the same as the sum of magnetic and viscous dissipation in 
simulations with nonzero v and rj. The runs 256H-£? 1 and 144H-£? 1 are exactly the same 
as the runs 256H-£? 1 and REF2 in CV00. The energy spectra as a function of time for these 
runs can be found in that paper. 



2.2. Defining the Local Frame 

The GS95 model deals with strong MHD turbulence and should be distinguished from 
theories that deal with weak MHD turbulence (e.g. Sridhar & Goldreich 1994; Ng & Bhat- 
tacharjee 1996; Galtier et al. 2000; see also Goldreich & Sridhar 1997). In strong MHD 
turbulence eddy-like motions mix up magnetic field lines perpendicular to the local direction 
of magnetic field. Thus, as in the case of hydrodynamic turbulence, the correlation time for 
coherent structures is comparable to the inverse of k±Vk for any scale kj 1 . These mixing 
motions are strongly coupled to wave-like motions with a correlation time (/c||Va) _1 . The 
GS95 model is based on the concept of a 'critical balance' between these time scales, that 
is k\\VA ~ k±Vk- ■ This results in a scale dependent anisotropy, h\\ oc k 2 / 3 so that the eddies 
are increasingly elongated on smaller scales. 

The turbulent magnetic field changes its direction in the global system of reference. It is 
important that the mixing motions are available only in the direction perpendicular to the 
local direction of magnetic field. Thus the theory must be formulated using the system of 
reference aligned with the local magnetic field. CV01 discusses in detail one way of defining 
this system given numerical data. 
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Figure 1 is a schematic representation of the GS95 model. 1 In Fourier space, the energy 
injected on large scales excites large scale Fourier components of the magnetic field (the 
dark region at the center in Fig. la). The external magnetic field makes the subsequent 
energy cascade to small scales anisotropic: it occurs in the directions perpendicular to the 
mean external field. The GS95 model states that most of the energy is confined to the 
region k\\ = ±k 2 / 3 , and as the energy cascades to larger values of k± the energy of Fourier 
components between k± and k± + 1 decreases as E(k±) oc k^^ . 

As illustrated in Fig. lb, eddies are not aligned along the mean field Bo- Instead, 
they are aligned along the local mean field lines. The local mean magnetic field defines the 
physically relevant background for the eddy dynamics, and is determined by the Fourier 
components whose wavenumbers are a bit less than the characteristic wavenumber of the 
eddy. In practice, it can be obtained by averaging the magnetic field in the vicinity of 
the eddy over a volume slightly (say, for example, 2x) larger than the size of the eddy 
(see CVOO for details). The solid curves in Fig. lb represent this kind of locally defined 
mean, formed by all magnetic Fourier components whose scales are a bit larger than eddy 
1 (or 1'). The characteristic scale of this wandering is L ~ 1/fcz,, the energy injection scale, 
because eddy 1 (or 1') is only slightly smaller than the energy injection scale. This large- 
scale wandering is smooth but dominates over smaller scale effects because the the magnetic 
energy is concentrated on larger scales. Wandering by smaller scale magnetic fields is weaker 
and causes smaller deviations from the large-scale wandering. We depict the additional 
wandering caused by scales a bit larger than the eddy 2 as a dashed curve in Figure lb. For 
the eddy 1, the solid curve defines the local mean magnetic field and, for the eddy 2 the 
dashed curve. 

The GS95 model is dominated by local dynamics, that is, in this model disturbances 
lose their coherence when propagating over a single wavelength. To the extent that the 
dynamics are local it is obvious that the only relevant magnetic field is the local mean field. 
As an example, consider eddies 1 and 2 in Figure lb again. For eddy 1, the solid curve can 
be regarded as a local mean field line and the energy cascade will take place perpendicular 
to this field, along the direction AA' . Smaller eddies, like eddy 2 sees the dashed curve as 
its local background magnetic field with a slightly different direction for the energy cascade. 
Since the difference between the two curves is small, the direction of the energy cascade 
differs only slightly as a function of scale. Basically, energy cascades along AA' in the region 
near the eddy 1. Similarly, energy cascades along DD' in the region near the eddy 1'. We 
are left with an energy cascade which differs both as a function of scale and as a function 



Note that Figure 1 is local frame representation of the GS95 model. See below for more details. 
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Fig. 1. — (a) Fourier space structure, (b) Real space structure. Large eddies (eddy 1 or 
1')) have similar semi-major axis (~ I/Zcim) and semi-minor axis (~ l/^J-,i)- Therefore, 
they are almost isotropic. Smaller eddies (eddy 2 or 2') have relatively larger semi-major 
axis (~ 1/^11,2) to semi-minor axis (~ l/fcj,^) ratio. Therefore, they are relatively more 
elongated. Energy cascades in the directions perpendicular to large-scale magnetic field 
lines (e.g. direction AA' or DD' in both figure (a) and (b)). This effect obscures scale- 
dependent anisotropy (e.g. fcn oc k 2 / 3 in GS95 model) when we perform Fourier analysis in 
the global frame. In Figure (b), the solid curves represent the wandering of magnetic field 
lines by the large scale magnetic fields. The solid curves can define the directions of the 
local mean magnetic field line for eddy 1 or 1'. Similarly, the dashed curves can define the 
directions of the local mean magnetic field line for the eddy 2 or 2'. 
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of location. However, the dynamics are not entirely local, in the sense that disturbances do 
propagate along field lines without retaining phase coherence. Consequently, perpendicular 
motions in eddy 1 can lead to similar motions in eddy 1', even though AA' and DD' are not 
parallel vectors. This is one of the key features of the GS95 model, and carries with it the 
implication is that dynamic variables need to be evaluated in terms of the local mean field 
direction, rather than a global coordinate system. Conversely, the ability to generate a more 
meaningful description in terms of a local, and scale-dependent, coordinate system can be 
taken as an indirect confirmation of some features of the GS95 model. 

In summary, when we want to describe the scaling of eddy shapes, we should correctly 
identify the direction of the local mean magnetic field. When we talk about anisotropy, we 
talk about anisotropy with respect to local mean magnetic field lines. Because of this, it is 
necessary to introduce a 'local' frame in which the direction of the local mean magnetic field 
lines is taken as the parallel direction. When we consider the GS95 picture (i.e. k\\ oc k 2 / 3 ) 
in Fourier space, we are considering the local frame in real space, and vice versa. When 
we describe turbulence with respect to the global frame, which is fixed in real space, the 
corresponding Fourier space structure no longer shows the GS95 picture. Instead, we have 
a relation close to ku oc k±. This is because when energy cascades along AA', DD', or 
some intermediate directions in real space (Fig lb), it cascades along the directions between 
AA' and DD' in Fourier space (Fig la), which implies that, when we perform the Fourier 
transform with respect to the fixed global frame 2 , we will get k\\ oc k±. The true scaling 
relation is eclipsed by the wandering of large scale magnetic field lines. 

It is very important to identify the local frame. In this paper, when we calculate decay 
time scale, intermittency, and the correlation tensor, we always refer to the local frame. 

3. Scaling Relations 

3.1. Time Scale of Motions 

One of the basic questions in the theory of MHD turbulence is the slope of the one- 
dimensional energy spectra. As we have seen, GS95 obtained a spectral index of —5/3. 
In the numerical simulations of CVOO the spectral index is close to —5/3, while it is very 
close to —3/2 in MG01. The IK theory predicts a k~ 3 ^ 2 scaling, although the other features 
of this model are definitely inconsistent with all the numerical evidence. MG01 attributed 



2 Many astrophysical observations, for instance, interferometric observations of turbulent HI (Lazarian 
1995), provide the statistics measured in the global frame. 
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their result to the appearance of strong intermittency in their simulations. We note that 
the inertial range of the solar wind shows a spectral index of —1.7 (Leamon et al. 1998; see 
also Matthaeus & Goldstein 1982), but this number should be considered cautiously. The 
physics of the solar wind is undoubtedly more complicated than the simulations described 
here. 

Can we test which scaling is correct? The cascade time as function of scale presents us 
with an interesting constraint. 

IK theory and GS95 model predict different scalings for turbulent cascade time scale 
(t cas ). In both theories, t cas can be determined by the scale-independence of the cascade: 

v l/tcas = constant. (6) 

Since v\ is proportional to kE(k), we have 

tcas,IK OC k~ 1/2 , tcaa,GS OC k~ 2/3 (7) 

for IK theory and GS95 model respectively. This result is also useful for certain intermittency 
theories (see §3.3). MG01 studied the cascade time scale using three different methods and 
obtained slopes comparable to —2/3 (i.e. t cas oc k ± 2 ^ 3 ) in two methods and —1/2 in the 
other method. 

Here we consider a different method of evaluating t cas . The purpose of our calculation is 
to test MGOl's result using another numerical method and demonstrate the effects of large 
scale fluid motions on the calculation of t cas . 

Symbolically, we can rewrite the MHD equations as follows: 

v k = iVZ, (8) 
bk = Nl (9) 

where N v and N b represent the nonlinear terms. We have ignored the dissipation terms. 
Naively speaking, we might obtain the time scale by dividing |v k | by \N£\. However, this 
gives t cas oc fc -1 , where the exponent is almost exactly —1. This is not actually a measure 
of the cascade time. We note that CVOO obtained a similarly misleading relation for the 
cascade time, and attributed it to the effect of large scale translational motions. Although 
they used a different method to calculate the cascade time, the same argument applies here. 
If we consider the interaction between a small eddy and a large scale (translational) fluid 
motion, then the translation can be removed by a Galilean transformation, and there is no 
associated energy cascade. However, the phase of the Fourier components that represent 
the small eddy is affected by the large scale translational motion, and changes at a rate 
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kV, where V is the large scale velocity. The corresponding nonlinear term has a magnitude 
of |iVfc| ~ IvfclfcV, which accounts for the (misleading) relation t cas ~ fc -1 . The cascade 
time as a function of wavenumber can be evaluated directly from our simulations, but only 
after we filter out translational motions arising from eddies much larger than the scale under 
consideration. 

We correct for the presence of large scale motions by restricting the evaluation of the 
nonlinear terms to contributions coming from the interactions between the mode at k and 
other modes within the range 3 oik/2 and 2k. In doing this, we retain the uniform magnetic 
component Bq. We show the result in Figure 2. Our result supports the GS95 model: 
t cas oc A; -2 / 3 . In comparison with MG01 we obtained this result using a different method and 
for a different kinetic/magnetic energy ratio. 

In the GS95 model, t cas is determined by the relation t cas ~ l±/vi±. This means that 
the cascade time scale is virtually synonymous with the eddy turnover time, which is also 
true for hydrodynamic turbulence. It is obvious that the cascade time determines the decay 
time scale of turbulence. As a consequence, the GS95 model implies that MHD turbulence 
decays as fast as hydrodynamic turbulence (say, t decay = a few eddy turnover times). Note 
that no matter how strong the external field is, strong MHD turbulence decays within a few 
eddy turn over times. We discuss the implications of this result, and some limitations, in §5. 

These results support the original GS95 theory. However, we are not in a position to 
directly confront the results of MG01. Our simulations differ from theirs in many ways, 
including the shape of the computational box, the range of length scales, and the strength 
of the uniform background field. We shall address those issues elsewhere. In §3.3, we will 
discuss what the study of intermittency implies about the slope of energy spectra. 

3.2. Velocity Scaling 

In the GS95 model, v\ is proportional to ~ k±E(k±), where E(k±) is the one-dimensional 
energy spectrum. Since E(k±) oc k^ 3 , we have v \ oc I 2 / 3 or SF 2 (l±, 0) oc I 2 / 3 , where SF 2 is 
the second order structure function: 

^ 2 (/±,/||)=<[v(x + l)-v(x)] 2 >, (10) 

where the angle brackets denote the spatial average over x and 1 = /j_e_i_ + l\\e^. The vectors 
ej_ and ey are unit vectors perpendicular and parallel to B £ respectively. The vector B L 



3 This assumes some sort of locality which may not be exact in the presence of strong intermittency. 




Fig. 2. — Cascade time scale. The GS95 model predicts t cas oc k 2 / 3 , while the IK theory 
predicts t cas oc k~ l l 2 . Our result supports the GS95 model. Run 256H-_B 1- 

Fig. 3. — Second-order structure functions. Across local mean magnetic field lines, the 
second-order structure functions follow r 2 / 3 . Along the local mean magnetic field lines, they 
follow r 1 . For pseudo-Alfven modes this defines the scaling of motions parallel to local mean 
magnetic field lines. V SF 2 and BSF 2 denote the second order structure functions for velocity 
and magnetic field respectively. Run 256H-fi l- 
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denotes the local mean magnetic field. 

What can we say about the velocity scaling parallel to B L ? We can consider two 
different quantities. First, we can consider the scaling of Alfven components in the direction 
parallel to B L . Second, we can also consider the scaling of pseudo-Alfvenic components 
along the direction of B^. The pseudo-Alfvenic components are incompressible limit of slow 
magnetosonic waves. While they have the same dispersion relations as the Alfven modes, 
the velocities Vi of polarization are completely different: the directions lie in the plane 
determined by B and k. Note that the Alfven modes have velocities perpendicular to the 
plane determined by B and k. 

There are several ways to derive the scaling relation for Alfvenic turbulence from the 
GS95 model. First, suppose that the second order structure function along local B^ follows 
a power law: ST^O,/]]) oc Z™. When we equate SF 2 (0,l\\) and SF 2 (lj_, 0), we should retrieve 

the GS95 scaling relation, oc I 2 / 3 (see MG01). We conclude that m = 1 and SF 2 (0, oc /y. 
Since k\\E v {k\\) (oc v 2 ) oc SF 2 (l\\,0), we have E v {k\\) oc k^ 2 . Alternatively, we can write 

El{k u k\\)^kl l0/3 g{k\\/kf ), (11) 

where E^(k±, k\\) is the 3-dimensional energy spectrum and g is a function which describes 
distribution of energy along k\\ direction in Fourier space. We will give a reasonable fit to its 
functional form in the next section. 

We plot our results in Figure 3, in which we observe that SF 2 (l±, 0) (across B^) oc I 2 / 3 
and 57*2(0, Z||) (parallel to B^) oc The velocity field follows these relations quite well, while 
the magnetic field follows a slightly different relation across B^. Both Alfvenic and pseudo- 
Alfvenic components follow similar scalings in the directions parallel to B^ In section 4 
we shall also show that the 3D spectrum of pseudo-Alfven motions has a form similar to 
(11). On this basis we conclude that the scaling E v {k\\) oc k^ 2 also applies to pseudo-Alfven 
motions, where the velocities are mostly parallel to the local mean magnetic field B L . The 
corresponding velocities v 2 oc k\\E v {k\\) which means that v 2 oc k^ 1 oc This result is 
important for many problems, including dust transport (Yan et al 2001). Note that the 
energy spectrum is steeper when expressed as a function of the parallel direction. 

In this subsection we extended the GS95 model for the parallel motions and pseudo- 
Alfven modes and confirmed it through numerical simulations. 
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3.3. Intermittency 

MG01 studied the intermittency of dissipation structures in MHD turbulence using the 
fourth order moments of the Elsasser fields and the gradients of the fields. Their simulations 
show strong intermittent structures. We use a different, but complementary, method to 
study intermittency, based on the higher order longitudinal structure functions. Our result 
is that by this measure the intermittency of velocity field in MHD turbulence across local 
magnetic field lines is as strong, but not stronger, than in hydrodynamic turbulence. 

In fully developed hydrodynamic turbulence, the (longitudinal) velocity structure func- 
tions S p =< ([v(x + r) — v(x)] -r) p >=< 5v p L (r) > are expected to scale as r** p . For example, 
the classical Kolmogorov phenomenology (K41) predicts ( p = p/3. The (exact) result for 
p=3 is the well-known 4/5-relation: < 8v p L (r) >= — (4/5)er, where e is the energy injection 
rate (or, energy dissipation rate). On the other hand, She & Leveque (1994, hereinafter S-L) 
proposed a different scaling relation: (p L = p/9 + 2[1 — (2/3) p / 3 ]. Note that She-Leveque 
model also implies ( 3 = 1. 

So far in MHD turbulence, to the best of our best knowledge, there is no rigorous 
intermittency theory which takes into account scale-dependent anisotropy. Therefore, we 
will use an intermittency model based on an extension of a hydrodynamic model. Politano 
& Pouquet (1995) have developed an MHD version of She-Leveque model: 

C p = - g {i - x) + C (1 - (1 - x /cy'») , (12) 

where C is the co-dimension of the dissipative structure, g is related to the scaling 
and x can be interpreted as the exponent of the cascade time t cas oc l x . (In fact, g is related 
to the scaling of Elsasser variable z: z\ ~ l 1 ^ 9 .) In the framework of the IK theory, where 
g — 4, x — 1/2, and C = 1 when the dissipation structures are sheet-like, their model of 
intermittency becomes (IK =p/8 + l- (l/2) p / 4 . On the other hand, Miiller & Biskamp 
(2000) performed numerical simulations on decaying isotropic MHD turbulence and obtained 
Kolmogorov-like scaling (E(k) ~ A;~ 5 / 3 and t ~ Z 2 / 3 ) and sheet-like dissipation structures, 
which implies g — 3, C — 1, and x = 2/3. From equation (12), they proposed that 

C B =p/9 + l-(l/3) p / 3 . (13) 

How does anisotropy change intermittency? We have determined the scaling exponents 
numerically, working in the local frame. We performed a simulation with a grid of 144 3 
collocation points and integrated the MHD equations from t=75 to 120. We calculated the 
higher order velocity structure functions for 75 evenly spaced snapshots. We average over 5 
consecutive values since the correlation time of the turbulence corresponds to 5 snapshots. 
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We calculated the scaling exponents from these averaged structure functions. We obtained a 
total of 15 (=75/5) such structure functions and scaling exponents. We believe these 15 data 
sets are mutually independent. We plot the result in Figure 4a. The filled circles represent the 
scaling exponents of longitudinal velocity structure functions in directions perpendicular to 
the local mean magnetic field. It is surprising that the scaling exponents are close the original 
(i.e. hydrodynamic) S-L model. This raises an interesting question. In our simulations, we 
clearly observe that t cas oc I 2 / 3 and E(k) oc A;~ 5 / 3 . It is evident that MHD turbulence 
has sheet-like dissipation structures (Politano, Pouquet & Sulem 1995). Therefore, the 
parameters for our simulations should be the same as those of Miiller & Biskamp's (i.e. g = 3, 
C — 1, and x = 2/3) rather than suggesting C — 2. We believe that this difference stems 
from the different simulation settings: their turbulence is isotropic and ours is anisotropic. 
In fact, we expect the the small scale behavior of MHD turbulence should not depend on 
whether or not the largest scale fields are uniform or have the same scale of organization as 
the largest turbulent eddies. Nevertheless, given the limited dynamical range available in 
these simulations, it would not be surprising if the scale of the magnetic field has a dramatic 
impact on the intermittency statistics. It is not clear how scale-dependent anisotropy changes 
the intermittency model in equation (12) and we will not discuss this issue further. Instead, 
we simply stress that we have found a striking similarity between ordinary hydrodynamic 
turbulence and MHD turbulence in perpendicular directions. MG01 attributes the deviation 
of their spectrum from the Kolmogorov-type to the turbulence intermittency present in the 
MHD case. Since we do not reproduce their power spectrum, the fact that our intermittency 
statistics do not support this conjecture is unsurprising. Clearly more studies of the issue 
are necessary. 

In figure 4a, we also plot the scaling exponents (represented by filled squares) of lon- 
gitudinal velocity structure functions along directions of the local mean magnetic field. 
Although we show only the exponents of longitudinal structure functions, those of trans- 
verse structure functions follow a similar scaling law. It is evident that intermittency 
along the local mean magnetic field directions is completely different from that of previ- 
ous (isotropic) models. Roughly speaking, the scaling exponents along the directions of 
local magnetic field are 1.5 times larger than those of perpendicular directions. Interestingly 
this result implies anisotropy becomes scale independent under the following transforma- 
tion: (r±,r\\) — > (r±,r?/ 3 ). This is consistent with the idea that eddies are stretched along 
the directions of the local mean magnetic field - if we shrink them in the scale-dependent 
manner described above along the local field lines, the result is similar to ordinary hydro- 
dynamic turbulence. In this interpretation it is not surprising that MHD turbulence looks 
similar to ordinary hydrodynamic turbulence across the local mean magnetic field lines - the 
scaling relation in perpendicular directions is not affected by the local mean magnetic field. 
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Clearly this result is hard to explain using previous models, for example the IK theory. The 
error bars are larger for parallel directions because fewer number of pairs are available for 
calculation of the structure functions in these directions than perpendicular directions. 

In Figure 4b, we plot average longitudinal velocity structure functions. The slope of 
the third order structure function is very close to 1. The third order structure function is 
slightly different from the one discussed earlier in that we calculate < |5fx| 3 (r) > instead of 
< (^L) 3 (r) >. 

The second order exponent ( 2 is related to the the 1-D energy spectra: E(k±) oc k^ 1+< * 2 \ 
Previous 2-D driven MHD calculations for Bq = by Politano, Pouquet, & Carbone (1998) 
also found £2 ~ 0.7. However, Biskamp, & Schwarz (2001) obtained (2 ~ 0.5 from decaying 
2-D MHD calculations with B = 0. Our result suggests that (2 is closer to 2/3, rather 
than to 1/2. (It is not clear whether or not the scaling exponents follow the original S-L 
model exactly. At the same time, our calculation shows that the original S-L model can be 
a good approximation for our scaling exponents. The S-L model predicts that (2 ~ 0.696.) 
Therefore, our result supports the scaling law E(k±) oc kj 5 ^ 3 at least for velocity. For the 
parallel directions, the results support E(k\\) oc k7 2 although the uncertainty is large. 



4. The MHD Fluctuation Tensor 



For many purposes, e.g. cosmic ray propagation and acceleration, heat transfer etc., it 
is necessary to know the tensor describing the statistics of magnetic and velocity field. For 
those applications the one-dimensional spectrum described in MG01 is not adequate and a 
more detailed description is necessary. 

General second-rank correlation tensors are important tools in the statistical description 
of turbulence. Oughton, Radler, & Matthaeus (1997) gave a comprehensive formalism for 
the tensors for MHD turbulence and we use their results as a starting point of our argument. 
Consider the velocity correlation tensor 



--<v i (x)v j (x. + r)>, (14) 

where the angle brackets denotes an appropriate ensemble average. The Fourier transform 
of this tensor is 

=< «,(k)€}(k) >, (15) 

where the asterisk denotes complex conjugate. We can rewrite equation (20) of Oughton et 
al. (1997) as 



QV 



kikj 



E" k + 



k 2 



ek) 



F v (k)+X ij: (16) 
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where E v (k) is (3D) kinetic energy spectrum of all (shear + pseudo) modes, F v (k) is the 
difference of shear- Alfven energy and pseudo-Alfven energy at wave vector k divided by k"j_, e 
is a unit vector along B , and = —i[5i^€j a p + 5j^e ial 3}e a kp(e^k 2 — k^e-]i)C v +i€ij a k a H v is a 
term that describes deviation from mirror symmetry. In this paper we consider axisymmetric 
turbulence (caused by Bo) with mirror symmetry so that Xij = 0. We need only two 
scalar generating functions, E v and F v , for the correlation tensor. This is consistent with 
Chandrasekhar (1951; see also Oughton et al. 1997). 

In this subsection, we will show that the tensor is suitably described by 



Sij — Ai 



k{kj 

5 - - -w 



k 



-10/3 



_L 



exp — A 2 



1.2/3 



(17) 



where A 1 ~ Bq/L 1 ^ 3 and A 2 ~ L 1 ^ 3 are parameters. 

First, we choose e=(0,0,l), the direction of B . Then, equation (16) becomes 



QV 



E v (k) + 



kk- 

[5 3 ikj + 5 3j ki]k 3 - 5 3i 5 3j k 2 l -rJ-k 2 



k 



2 ^3 



F v (k). 



(18) 



In the absence of anomalous damping of the pseudo-Alfven modes, as in our simulations, 
we can show that F in above expression is negligibly small. Note that F = (S — P)/k\, 
where S and P are the squares of the amplitudes of the shear- and pseudo-Alfven modes (i.e. 
3D energy spectra). To evaluate S and P, we measured their strength in global frame. (It 
is non-trivial to correctly define Alfven modes and pseudo-Alfven modes in the local frame.) 
Figure 5 shows that they have similar strengths. We assume that the same relation holds 
true in the local frame. Since F is the difference between S and P, it follows that F(k) is 
small compared with E(k). 

In the previous paragraph, we assumed that there is no special damping mechanism 
for the pseudo-Alfven modes. However, it is known that pseudo-Alfven modes in the ISM 
are subject to strong damping due to free streaming of collisionless particles along the field 
lines (Barnes 1966; Minter & Spangler 1997). When the pseudo-Alfven modes are absent, 
equation (16) becomes 



k{kj 



E v (k) + 



[S 3i kj + S 3j ki]k 3 - S 3i S 3j k 2 - 



kjkj 2 

~W 3 



E v (k) 



(19) 



For i, j = 1, 2, this becomes = (5ij — kikj/k 2 )E — (kikj/k 2 )(k\\/k±) 2 E fx (5ij — kikj/k 2 )E. 
And, it is easy to show that S i3 = S 3i = 0. This is easily understood when we note that 
shear- Alfven waves do not have fluctuations along B . 
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In summary, the tensor reduces to 

(1 - k\jk 2 )E -{hk 2 /k 2 )E -{hh/k^E 

-(hk 2 /k 2 )E (1 - k\fk 2 )E -(k 2 k 3 /k 2 )E I (20) 

-(hh/k^E -(k 2 k 3 /k 2 )E (1 - k\jk 2 )E 

for turbulence with both Alfvenic and pseudo-Alfvenic components, and 

(1 - k\jk 2 )E -(k x k 2 /k 2 )E 

-(hk 2 /k 2 )E (1 - k\jk 2 )E | (21) 


for shear Alfvenic turbulence. In equation (21), E stands for the energy of Alfven components 
only, which is roughly one half of the E in equation (20). 

The remaining issue is the form of E. Note that the trace of S^j is 2E V . In real space, 
the trace is the velocity correlation function. Consequently, we can obtain E v through a 
FFT of the real-space velocity correlation function, which is directly available from our 
data cube. However, the velocity correlation function in real space contains considerable 
numerical noise. In order to minimize its effects while obtaining an empirically useful form 
for E we first guess E v in Fourier space, do the FFT transform, and then compare the 
transformed result with the actual velocity correlation function. Since the trace of is the 
(3-dimensional) energy spectrum in Fourier space, we start with the original expression in 
Goldreich & Sridhar (1995) given by equation (11): 

E ^h)~^-^p, (22) 

where the functional form of g(y) was not specified. We have tried several functional forms 
for g ; Gaussian (g oc exp(— Bk 2 /k 4 / 3 )), exponential (exp(— Bk\\/k 2 / 3 )), and a step function. 
We have found that an exponential form for g gives the best result (Figure 6). Fig 6a is the 
actual data we obtained from our simulation and Fig 6b is the Fourier-transformed velocity 
correlation function. Note the similarity of the contours in both plots. We conclude that the 
tensor can be suitably described by equation (20) or (21) with 

E(k ± , k\\) = (B /L^)k- ± w/3 exp , (23) 



where k± and k\\ should both be interpreted as the absolute magnitudes of those wavevector 
components. 
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However, it is worth noting a clear limitation of equation (23): it has a discontinuous 
derivative near k\\ = 0. One way to overcome this difficulty is to use the Castaing function 
(Castaing, Gagne & Hopfinger 1990) 4 

which is smooth near zero but looks exponential over a broad range. It is possible to see 
that for A = I and a = k 2 / 3 /L^ 3 , 

( n/z k \\\ 1 f°° ( k \\\ ( \n\L^a/kT)\da 
exp [-L^ y-J Q exp j exp \- j - (25) 

However, for many practical applications, we feel that the expression in equation (23) is 
adequate. For instance, in a forthcoming paper (Yan et al. 2001) this tensor is used for 
describing cosmic ray propagation and we find a strong suppression of cosmic ray scattering 
compared with the generally accepted estimates (for example, Schlickeiser 1994). However, 
if the behavior around k\\ is important, the Castaing function would be preferred. In our 
simulations there is no way to distinguish between exponential and Castaing distributions. 



5. Decay of MHD Turbulence 

Turbulence plays a critical role in molecular cloud support and star formation and the 
issue of the time scale of turbulent decay is vital for understanding these processes. 

If MHD turbulence decays quickly then serious problems face researchers attempting to 
explain important observational facts, i.e. turbulent motions seen within molecular clouds 
without star formation (see Myers 1999) and rates of star formation (Mckee 2000). Earlier 
studies attributed the rapid decay of turbulence to compressibility effects (Mac Low 1999). 
Our present study, as well as earlier ones (CV00, MG01), shows that turbulence decays 
rapidly even in the incompressible limit. This can be understood in the framework of GS95 
model where mixing motions perpendicular to magnetic field lines form hydrodynamic-type 
eddies. Such eddies, as in hydrodynamic turbulence, decay in one eddy turnover time. 

How grave is this problem? Some possibilities for reconciling theory with observations 
were studied earlier. For instance, some problems may be alleviated if the injection of energy 



4 Our motivation for introducing the Castaing function is phenomenological, not theoretical. The theory 
of the distribution function for MHD turbulence is uncertain, and far beyond the scope of this paper. 
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happens on the large scale, the eddies are huge, and the corresponding time scales, are much 
longer (see Lazarian 1999). The fact that the turbulence decays according to a power-law, 
rather than exponentially, also helps. Indeed, if turbulent energy decays as t -1 , as suggested 
by Mac Low wt al. (1998), a substantial level of turbulence should persist after 4-5 turnover 
times. 

There is, however, another property of astrophysical turbulence related to the peculiar 
nature of the energy injection. It is accepted that sources of interstellar turbulence are 
localized. As a result, there is a substantial imbalance between the ingoing and outgoing 
energy flux surrounding every source. Below we consider the effect of this imbalance on the 
turbulence decay time scale. 

For an imbalanced turbulence, it is useful to consider the Elsasser variables, = v±b, 
which describe wave packets traveling in opposite directions along the magnetic field lines. 
Imbalanced turbulence means that wave packets traveling in one direction (say, z + ) have sig- 
nificantly larger amplitudes than the other. In astronomy, many energy sources are localized. 
For example, SN explosions and OB winds are typical point energy sources. Furthermore, 
astrophysical jets from YSOs are believed to be highly collimated. With these localized 
energy sources, it is natural to think that interstellar turbulence is typically imbalanced. 
In fact, the concept of an imbalanced cascade is not new. Earlier papers (e.g., Matthaeus, 
Goldstein & Montgomery 1983; Ting, Matthaeus & Montgomery 1986; Ghosh, Matthaeus 
& Montgomery 1988) have addressed the role and evolution of cross-helicity (= v • b). Since 
4 < v ■ b >=< (z + ) 2 > — < (z~) 2 >, non-zero cross-helicity implies an imbalanced turbu- 
lent cascade. These works, however, were mainly concerned with the growth of imbalance 
in decaying turbulence. Ghosh et al. (1988) investigated the evolution of cross-helicity and 
various spectra in driven turbulence. Hossain et al. (1995) discussed the effects of cross 
helicity and energy difference D =< v 2 > — < b 2 > on the decay of turbulence. Their low 
resolution 3D numerical simulations show the effect of cross helicity, although that effect 
is not very conspicuous. A further study of imbalanced turbulence was given in MG01, 
who also suggested a connection between spontaneous appearance of local imbalance in the 
turbulent cascade and intermittency in MHD turbulence. 

In this subsection, we explicitly relate the degree of imbalance and the decay time scale 
of turbulence in the presence of a strong uniform background field. 

In Figure 7 we demonstrate that an imbalanced cascade does extend the lifetime of 
MHD turbulence. We use the run 144H-_B 1 to investigate the decay time scale. We ran the 
simulation up to t=75 with non-zero driving forces. Then at t=75, we turned off the driving 
forces and let the turbulence decay. At t=75, there is a slight imbalance between upward and 
downward moving components (E + = 0.499 and E_ = 0.40). This results from a natural 
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fluctuation in the simulation. The case of (E + ) to = 80%(E) tQ corresponds to the simulation 
that starts off from this initial imbalance. In other cases, we either increase or decrease the 
energy of z~ components and, by turning off the forcing terms, let the turbulence decay. We 
can clearly observe that imbalanced turbulence extends the decay time scale substantially. 
Note that we normalized the initial energy to 1. The y-axis is the total (=up + down) 
energy. 

In Figure 7b, we re-plot the Figure 7a in log-log scale. For balanced case (i.e. zero cross- 
helicity case; solid curve), the energy decay follows a power law E(t) oc t a , where a is very 
close to 1. This result is consistent with previous 3D result by Hossain et al. (1995). Note 
that hydrodynamic turbulence decays faster than this. For example, Kolmogorov turbulence 
decays as E(t) oc t~ 10 ^ 7 . In this sense, it may not be absolutely correct to say that both 
hydro and MHD turbulence decay within one eddy turnover time. However, note that the 
power law does not hold true from the beginning of decay. We believe that, at the initial 
stage of decay, the speed of decay is still roughly proportional to the large-scale eddy turn 
over rate. 

How far does a wave packet travel when there is an imbalance? Consider the equations 
governing an imbalanced cascade. From the MHD equations, Hossain et al. (1995) and MG01 
derived a simple dynamical model for imbalanced turbulence. For decaying turbulence, they 
found 

" E + - (26) 



dt 

dE_ E_E_ 



1/2 



it L ' < 27 > 

where L is the largest energy containing eddy scale. 5 From these coupled equations, they 
showed that imbalance grows exponentially in decaying turbulence. Now let us consider a 
large amplitude wave packet traveling in an already (weakly) turbulent background medium. 
Suppose that the large amplitude wave packet corresponds to z + . Using the simplified 
equations, we obtain E + /E + = —E x l 2 jL. If the background turbulence has a constant 
amplitude, the z + wave decays exponentially. It can travel 

v + At ~ (E¥ 2 /El /2 )L, (28) 

where we use v + ~ El /2 . This means that the wave packet can travel a long distance when 

1 /2 1 /2 

imbalance is large (i.e. E+ ^> E_ ). In real astrophysical situations, the problem is not as 



5 For our current purposes, these simple system of equations are enough. For more rigorous equations on 
the evolutions of E±(k), readers may refer to earlier closure equations, e.g. Grappin et al. (1982). See also 
Hossain et al. (1995) for time evolution of L. 
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simple as this. Instead, the wave packet and the background turbulence can have different 
length scales (as opposed to the single scale L in the equations). We also need to consider the 
fact that the amplitude of the background turbulence does not stay constant and the front 
of the wave packet decays faster than the tail of the packet. Finally, MHD turbulence can 
influence the pressure support only if turbulent motions are at least comparable to the sound 
speed, which obviously requires a fully compressible treatment. Preliminary calculations 
with a compressible code (Cho & Lazarian 2002) show marginal coupling of compressible 
and incompressible motions, while the development of the parametric instability (Fukuda & 
Tomoyuki 2000) requires more time to develop. We plan to investigate these possibilities in 
future. 

In this section, we found that turbulence decay time can be slow. This finding is very 
important for many astrophysical problems. 

6. Discussion 

How relevant are our calculations for the "big picture"? First of all, they provide more 
support for the GS95 theory, indicating that for the first time we have an adequate, if 
approximate, theory of MHD turbulence. Second, they extend the theory by treating new 
cases, e.g. an imbalanced cascade. Third, they establish new scaling relations and determine 
critical parameters, e.g. the functional form of g in equation (11), that will allow the theory 
to be applied to various astrophysical circumstances. 

Our calculations are made within an intentionally simplified model, which is based on 
the physics of an incompressible fluid. This surely raises the question of the applicability 
of our scaling relations and conclusions to realistic circumstances. There are situations 
where our scalings should be applicable. For instance, turbulence at very small scales is 
small-amplitude and therefore essentially incompressible. Processes that depend on the fine 
structure of turbulence, like scintillation, reconnection, and the propagation of cosmic rays 
of moderate energies should be well described using our results. 

If we consider the interstellar medium at larger scales, it is definitely compressible, 
has a whole range of energy injection/dissipation scales (see Scalo 1987), and the relative 
role of vortical versus compressible motions being unclear. Nevertheless, we believe that 
our simplified treatment may still elucidate some of the basic processes. To what extend 
this claim can be justified will be clear when we compare compressible and incompressible 
results. However, if we accept that fast and slow MHD modes are subjected to fast collision- 
less damping (see Minter & Spangler 1997) the remaining modes are incompressible Alfven 
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modes. Those should be well described by our model when turbulence is supersonic but 
sub-Alfvenic. Our preliminary results (Cho & Lazarian 2002) show that the coupling of the 
modes is marginal even in compressible regime. Incidentally, recent studies of turbulence of 
HI in both our Galaxy and SMC (Lazarian 1999, Lazarian & Pogosyan 2000, Stanimirovic 
& Lazarian 2001) show the spectra of velocity and density consistent with the Kolmogorov 
scalings 6 . 

Our approach is complementary to MG01. They studied turbulence in the regime when 
the magnetic energy is substantially larger than the kinetic energy at the energy injection 
scale. Physically their regime reflects better the properties of turbulence on small scales, 
where the magnetic energy is indeed dominant. In our calculations the kinetic energy is equal 
to the magnetic energy at the energy injection scale and therefore they reflect, for instance, 
what is happening in the interstellar medium at the large scales. Our results show that even 
on those scales GS scaling is applicable. This suggests that the astrophysical turbulence 
may be well tested not only via scintillations, that reflect properties of the turbulence on the 
small scales, but with other techniques, e.g. synchrotron emission. 



7. Summary 

Our findings can be summarized as follows: 

The energy cascade time scale at a length scale I (~ l/k) is proportional to Z 2 / 3 (/c~ 2 / 3 ), 
which is consistent with the prediction of the GS95 model and numerical simulations by 
MG01 who used a different method to obtain this scaling. In this respect MHD turbulence 
is similar to hydrodynamic turbulence. This scaling is distinctly different from the prediction 
of Iroshnikov-Kraichnan theory, t cas oc 

We found that velocity fluctuations in the direction parallel to local magnetic field follow 
a similar scaling for both Alfvenic and pseudo-Alfvenic modes. We determined that parallel 
motions due to pseudo-Alfven perturbations obey the following scaling: v\\ ~ k^ 2 . This 
finding is important for practical applications, e.g. for description of dust grain motion. 

To study intermittency we calculated higher order longitudinal velocity structure func- 
tions in directions perpendicular to the local mean magnetic field and found that the scaling 
exponents are close to Cp L = p/9 + 2[1 — (2/3) p / 3 ]. As this coincides with the She-Leveque 
model of intermittency in hydrodynamic flow we speculate that there may be more sim- 
ilarities between magnetized and unmagnetized turbulent flows than has been previously 



6 If density acts as a passive scalar its spectrum mimics that of velocity over the inertial range. 
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anticipated. 

We obtained correlation tensors which provide a good fit for our numerical results. 
These tensors are valuable for theoretical applications, e.g. to describe cosmic ray transport. 

We found that the rate at which MHD turbulence decays depends on the degree of 
energy imbalance between wave packets traveling in opposite directions. A substantial degree 
of imbalance can substantially extend the decay time scale of the MHD turbulence and the 
distance the turbulence can propagate from the source. 

The authors are thankful to Peter Goldreich for attracting our attention to the case 
of imbalanced cascade and for many valuable discussions. We also thank John Mathis for 
many useful suggestions and discussions. AL and JC acknowledge the support of NSF 
Grant...., and EV acknowledges the support of NSF Grant AST-0098615. This work was 
partially supported by National Computational Science Alliance under CTS980010N and 
AST000010N and utilized the NCSA SGI/CRAY Origin2000. 
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Fig. 4. — (a) Intermittency. Our result (denoted by the filled circles) suggests that MHD 
turbulence looks like ordinary hydrodynamic turbulence when viewed across the local field 
lines. SL represents the original She-Leveque model for ordinary hydrodynamic turbulence. 
IK and MB stand for the IK theory and the Miiller-Biskamp model respectively. Error bars 
are for 1-a level. Run 144H-B 1. (b) The second, third, and fifth order longitudinal velocity 
structure functions. These are structure functions averaged over the time interval (75,120). 
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Fig. 5. — Energy distribution of shear- and pseudo-Alfven waves. They have similar energy 
distributions in Fourier space (global frame). Run 256H-£> 1. Numbers by the contours are 
log 10 (£ 3 (k)). 

Fig. 6. — Upper panel: velocity correlation function from simulations. Lower panel: velocity 
correlation function generated using the tensor in equation (23). Parallel and perpendicular 
directions are taken with respect to local mean magnetic field. 
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Fig. 7. — (a) Decay of unbalanced turbulence. Imbalanced cascade can extend decay time. 
It is clear that the decay of turbulence depends on the degree of imbalance. Run 144H-5ol. 
(b) Same as (a). Log-log scale. Note that about 3 time units is one eddy turn over time. 
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Table 1. Parameters 



Run a 


N 3 


V 


V 


B 




144 3 


3.20 x 10- 28 


3.20 x 10- 28 


1 


256H-5 1 


256 3 


6.42 x 10" 32 


6.42 x 10- 32 


1 



a For 256 3 (or 144 3 ) grids, we use the notation 256X- 
Y (or 144X-Y), where X = H or P refers to hyper- or 
physical viscosity; Y = B 1 refers to the strength of 
the external magnetic fields. 



